/***************************************************************
                             zap2.c 
only the routines from the original zap.c that are different in
serial vs. parallel
***************************************************************/
#include <stdio.h>
#include <malloc.h>
#include <unistd.h>
#include <stdlib.h>
#include <assert.h>
#include "clam.h"

#include <netinet/in.h>
#include "timer.h"
#include "parallel.h"
#include "serial.h"

/*
 This will dump the matrix after parallel pair testing to check correctness 
#ifdef PARALLEL
#define DEBUG_PRINT_MATRIX
#endif
*/

extern int autosave_stats;

extern int Zrule_grand, cnt_sim_FN, cnt_sim_FP, unburry, cnt_nodes;
extern int cnt_cmp, total_cmp, cnt_overlap;
extern int Zbuild_matrix(), Greedy_make_Zset(int,int,int), Do_BestOf();
extern void quitCB(), cleanCB(), showClam2(), clearMsg(), init_CpM();
extern void unbury(), Start_timer(), Zmark_possible_buried(), OK_Zset();
extern void init_Zset(), End_timer(), Zproj_msg(), ZbuildCpM(), settime();
extern void big_wrap_up(), Zfree_zap(), Zsim_F_check(), Zscore_grand2();
extern int  ZscoreCpM(); 

static int cnt_cmp;
/* sness */
char stats_print1[256];
char stats_print2[256];
#define BUFSZ 120

/**************************************************************
                 Zscore_grand
to determine what clones have best matches and hence, start with in zcalc
Zrule_grand = 1 
**************************************************************/

void Zscore_grand(int i, int exponent)
{
   pthread_mutex_lock(&ZZ.matrix[i].lock);


   if (Zrule_grand)   ZZ.matrix[i].grand += (exponent * exponent);
   else         ZZ.matrix[i].grand = MaX(exponent, ZZ.matrix[i].grand);

   pthread_mutex_unlock(&ZZ.matrix[i].lock);

}

/****************************************************************
                       Def: Zbuild_contigs_zap
Complete build
****************************************************************/
void Zbuild_contigs_zap()
{
int i, c, save, cnt, j;
CLONE *clp;
char msg[100];
 int max=0, ctgcnt[4], ctg, sn_ctgcnt[10];
 char time_buf[BUFSZ];
 time_t ltime;
 struct tm *curr_time;

/* For sn_ctgcnt */
/* Bin       #'s */
/*   0         2 */
/*   1         3-5 */
/*   2      5-10 */
/*   3      10-20   */
/*   4      20-50 */
/*   5      50-100    */
/*   6      100-250    */
/*   7      250-500 */
/*   8      500-1000 */
/*  9      1000+ */

fpc_timer_t pairs, greed;
#ifdef DEBUG_PRINT_MATRIX
   FILE* debug_log;
#endif

   FILE* stats_log;
    char stats_pathname[256];
 
  save=cnt=0;
  ctgcnt[0]=ctgcnt[1]=ctgcnt[2]=ctgcnt[3]=0;

  for (j = 0; j<10; j++)
   sn_ctgcnt[j] = 0;

  if (graphActivate(g99)) quitCB();
  else cleanCB();
  if (contigs[0].count==0) {
     sprintf(ZBuf,"No singletons"); showClam2(); 
     return;
  }
  clearMsg();
  save = getcontigsmax();
  init_CpM(1); Outlog;
  strcpy(Fn,"Create");
  for (i= contigs[0].start; i != -1; i = clp->next)
  {
     clp = arrp(acedata, i, CLONE);
     if (clp->clone[0] == '!') continue;
     if (clp->match[0]!=' ') {
        unbury(i);
     }
  }
  
  if (!Zbuild_matrix(0)) return;
  
  initTimer(&pairs);
  startTimer(&pairs);
  if (parallel_node_creation() < 0)
  {
     stopTimer(&pairs);
     goto bail;
  }
  stopTimer(&pairs);

#ifdef DEBUG_PRINT_MATRIX
  if ((debug_log = fopen("matrix.output", "w")) == (FILE*)0)
  {
     perror("open matrix.output");
     exit(-1);
  }
  
  fprintf(debug_log, "Top matrix:\n");
  fprintf(debug_log, "  name:  %s\n", ZZ.name);
  fprintf(debug_log, "  size:  %i\n", ZZ.size);
  fprintf(debug_log, "  max:   %i\n", ZZ.max);
  fprintf(debug_log, "  num:   %i\n", ZZ.num);
  fprintf(debug_log, "  ctg:   %i\n", ZZ.ctg);
/*  fprintf(debug_log, "  type:  %i\n", ZZ.type); Appears to be unset anywhere */
  fprintf(debug_log, "  tol:   %i\n", ZZ.tol);
  fprintf(debug_log, "  band:  %i\n", ZZ.band);
  fprintf(debug_log, "  cut:   %g\n", ZZ.cut);
  fprintf(debug_log, "  cutOk: %g\n", ZZ.cutOK);
  fprintf(debug_log, "\n");

  for (i = 0; i < ZZ.size; i++)
  {
     int j;
     
     fprintf(debug_log, "Contig #%i:\n", i);
     fprintf(debug_log, "  set:        %i\n", ZZ.matrix[i].set);
/*     fprintf(debug_log, "  sparse:     %i\n", ZZ.matrix[i].sparse); Appears to be unset till zorder.c */
     fprintf(debug_log, "  cbQ:        %i\n", ZZ.matrix[i].cbQ);
     fprintf(debug_log, "  high_match: %i\n", ZZ.matrix[i].high_match);
     fprintf(debug_log, "  mtype:      %i\n", ZZ.matrix[i].mtype);
     fprintf(debug_log, "  pick:       %i\n", ZZ.matrix[i].pick);
/*     fprintf(debug_log, "  box:        %i\n", ZZ.matrix[i].box);   uncertain as to use */
/*     fprintf(debug_log, "  ebox:       %i\n", ZZ.matrix[i].ebox);  uncertain as to use */
     fprintf(debug_log, "  grand:      %i\n", ZZ.matrix[i].grand);
     fprintf(debug_log, "  save_grand: %i\n", ZZ.matrix[i].save_grand);
     fprintf(debug_log, "  leftb:      %i\n", ZZ.matrix[i].leftb);
     fprintf(debug_log, "  rightb:     %i\n", ZZ.matrix[i].rightb);
     /* uncertain what extra is */
     fprintf(debug_log, "  n_extra:    %i\n", ZZ.matrix[i].n_extra);
     fprintf(debug_log, "  max_extra:  %i\n", ZZ.matrix[i].max_extra);
     fprintf(debug_log, "  nbands:     %i\n", ZZ.matrix[i].nbands);
     fprintf(debug_log, "  n_cb:       %i\n", ZZ.matrix[i].n_cb);
     /* cb_index shouldn't be setup yet */
     fprintf(debug_log, "  cin:        %i\n", ZZ.matrix[i].cin);
   /* first varies based on insertion order - threading makes it random */
     
     for (j = ZZ.matrix[i].first; j != -1; j = ZZ.pool[j].next)
     {
        fprintf(debug_log, "   (%i, %i) = %i - %i\n",
           ZZ.pool[j].i1, ZZ.pool[j].i2,
           ZZ.pool[j].exponent, ZZ.pool[j].match);
     } 
  }
  
  fclose(debug_log);
#endif

  Zmark_possible_buried();

  printf("After marking buried\n");

  initTimer(&greed);
  startTimer(&greed);
  while(Greedy_make_Zset(ZSET, 0, 0) && Do_BestOf()) {
     ctg = getcontigsmax()+1;
     sprintf(ZZ.name, "Ctg%d", ctg);
     OK_Zset(0, ctg);
     sprintf(contigs[ctg].projmsg,"%-6.3f %3d",Zset[ZS].score, Zset[ZS].nolap);
     contigs[ctg].ctgQs = Zset[ZS].nolap;

     max = MaX(max, Zset[ZS].n_clone);
     if (Zset[ZS].n_clone > 25) ctgcnt[0]++;
     else if (Zset[ZS].n_clone > 10) ctgcnt[1]++;
     else if (Zset[ZS].n_clone > 2) ctgcnt[2]++;
     else ctgcnt[3]++;
    /*   sness */
    if (Zset[ZS].n_clone > 1000) sn_ctgcnt[0]++;
    else if (Zset[ZS].n_clone > 500) sn_ctgcnt[1]++;
    else if (Zset[ZS].n_clone > 250) sn_ctgcnt[2]++;
    else if (Zset[ZS].n_clone > 100) sn_ctgcnt[3]++;
    else if (Zset[ZS].n_clone > 50) sn_ctgcnt[4]++;
    else if (Zset[ZS].n_clone > 20) sn_ctgcnt[5]++;
    else if (Zset[ZS].n_clone > 10) sn_ctgcnt[6]++;
    else if (Zset[ZS].n_clone > 5) sn_ctgcnt[7]++;
    else if (Zset[ZS].n_clone >= 3) sn_ctgcnt[8]++;
    else sn_ctgcnt[9]++;

     cnt++;
           /* clears ALL values for ZS set, except for olap (Q) and score */
     for (c=0; c< Zset[ZS].n_clone; c++) {
         i = Zset[ZS].csort[c].i;
         if (ZZ.matrix[i].n_extra>0) ZFREE(ZZ.matrix[i].extra, "extra");
         if (ZZ.matrix[i].cb_index != NULL) ZFREE(ZZ.matrix[i].cb_index, "cb_index");
         ZZ.matrix[i].n_extra=0;
         ZZ.matrix[i].cb_index = NULL;
         arrp(acedata, ZZ.matrix[i].cin, CLONE)->oldctg = ctg; 
     }
     i = Zset[ZS].ctg;
     init_Zset(ZS, 1);
     Zset[ZS].ctg = i;
  }

stop: ;
  stopTimer(&greed);
  LastCB=NoCB;
  sprintf(ZBuf,"\n"); Outlog;
  big_wrap_up(0);

  /* sness */
  /* Save Statistics */
  if (autosave_stats || messQuery("Save statistics of Build?")) {
   if (strlen(dirName) > 0) 
     strcpy(stats_pathname,dirName);
   else 
     strcpy(stats_pathname,".");
   strcat(stats_pathname,"/");
   strcat(stats_pathname,fileName);
   strcat(stats_pathname,".stats");
    if ((stats_log =  fopen(stats_pathname, "a")) == (FILE*)0) {
     printf("ERROR: Opening stats file (%s)",stats_pathname);
    } 
   else {
     if (time(&ltime) != -1) {
      curr_time = localtime(&ltime);
      if (!strftime(time_buf, BUFSZ, "%R %a %d %h %Y", curr_time)) 
        fprintf(stderr, "Warning: Not enough room for time stamp.");
     }
     else time_buf[0] = '\0';

     /* Print to stats file */
     fprintf(stats_log,"---------------------------------------\n");
     fprintf(stats_log,"Total Clones = %d\n", arrayMax(acedata));
     fprintf(stats_log,"Date %s\n",time_buf);
     fprintf(stats_log,"Real time taken:\n");
     fprintf(stats_log,"  Pairs: %s\n", formatTimeval(getRealTime(&pairs)));
     fprintf(stats_log,"  Greed: %s\n", formatTimeval(getRealTime(&greed)));
     fprintf(stats_log,"FPC Version %s\n",VERSION);
     fprintf(stats_log,"%s",stats_print1);
     fprintf(stats_log,"%s",stats_print2);
     fprintf(stats_log,"Create %d contigs (from ctg%d to ctg%d)\n",
            cnt, save+1, save+cnt); 
     if (cnt==0) {
      fprintf(stats_log,"No contigs created.\n");
     }
     else {
      fprintf(stats_log,"Contig sizes\n");
      fprintf(stats_log,"singles  = %d\n",contigs[0].count);
      fprintf(stats_log,"doublets = %d\n",sn_ctgcnt[9]);
      fprintf(stats_log,"3-5      = %d\n",sn_ctgcnt[8]);
      fprintf(stats_log,"6-10     = %d\n",sn_ctgcnt[7]);
      fprintf(stats_log,"11-20    = %d\n",sn_ctgcnt[6]);
      fprintf(stats_log,"21-50    = %d\n",sn_ctgcnt[5]);
      fprintf(stats_log,"51-100   = %d\n",sn_ctgcnt[4]);
      fprintf(stats_log,"101-250  = %d\n",sn_ctgcnt[3]);
      fprintf(stats_log,"251-500  = %d\n",sn_ctgcnt[2]);
      fprintf(stats_log,"501-1000 = %d\n",sn_ctgcnt[1]);
      fprintf(stats_log,"1001+    = %d\n",sn_ctgcnt[0]);
     }
     fprintf(stats_log,"Markers = %d\n",arrayMax(markerdata));
     fclose(stats_log);
   }
  }

bail:
  if (cnt==0) {
     sprintf(ZBuf,"No contigs created.\n");
     Outlog;
  }
  else {
    sprintf(ZBuf,"Create %d contigs (from ctg%d to ctg%d)\n",
            cnt, save+1, save+cnt); Outlog;
    sprintf(ZBuf,
     "Contig sizes: Max %d, %d (>25), %d (25:10), %d (9:3), %d (=2), %d Singles\n",
      max, ctgcnt[0], ctgcnt[1], ctgcnt[2], ctgcnt[3], contigs[0].count); Outlog;
  }
  settime(&Bd.date);
  ZbuildCpM(); /* set table values for saving in fpc file */
  if (!Zbatch_flag) ClamCtgDisplay2();


  printf("Real time:\n");
  printf("  Pairs: %s\n", formatTimeval(getRealTime(&pairs)));
  printf("  Greed: %s\n", formatTimeval(getRealTime(&greed)));

}

/***************************************************************
                       Def: Zcreate_node
****************************************************************/
/** The pseudo-origonal Zcreate_node - not optimized worth a dime **/
/** For an optimized&parallelized&distributed full build see below **/
static int Zcreate_node(int i)
{
int  j, n;
ZNODE *node, *node2;
int exponent;
CLONE clonei, clonej;

     if (!Zbatch_flag && graphInterruptCalled()) {
          printf("F4 pre-mature termination of %s\n", Fn);
          return 0;
     }


    /* could be cemap yac with no bands */
     if (!loadCzfast(&C1z, ZZ.matrix[i].cin )) {
          ZZ.matrix[i].set  = ZNOT_USED;
          return 1;
     }
     clonei = arr(acedata, ZZ.matrix[i].cin, CLONE);
     cnt_nodes++;

     ZZ.matrix[i].cb_index = (int *) malloc((C1z.nbands+1)*sizeof(int));
     NOMEM3(ZZ.matrix[i].cb_index, "cb_index", 1);
     ZZ.matrix[i].cb_index[0] = ZONE;

     if (Zadd || Zorder) n= 0; else n = i+1;
     for (j = n; j < ZZ.size;  j++)
     {
         if (j==i) continue;
         if (ZZ.matrix[j].set==ZNOT_USED) continue;

         clonej = arr(acedata, ZZ.matrix[j].cin, CLONE);
         if (LastCB==Select && clonej.selected==0) continue;

         if (!loadCzfast(&C2z, ZZ.matrix[j].cin)) {
            ZZ.matrix[j].set  = ZNOT_USED;
            continue;
         }

         Zsulston(0);
         exponent = ZscoreCpM();


         if (exponent==0) { /* could be buried with poor overlap */
             if (ZZ.matrix[i].mtype == IAM_CHILD){
                if (strcmp(clonei.match, clonej.clone)==0) exponent=1;
             }
             else
             if (ZZ.matrix[j].mtype == IAM_CHILD){
                if (strcmp(clonej.match, clonei.clone)==0) exponent=1;
             }
         } 
         cnt_cmp++;
         total_cmp += 1.0;
         if (exponent == 0) continue;
         cnt_overlap++;

         Zscore_grand(i, exponent);
         node = Znew_node(i, j);
         node->match = Sz.match;
         node->exponent = exponent;

         Zscore_grand(j, exponent);
         node2 = Znew_node(j, i);
         node2->match = Sz.match;
         node2->exponent = exponent;

         if (clonei.mattype!=PSEUDO && strcmp(clonei.match, clonej.clone)==0) {
             ZZ.matrix[i].high_match = j;
             ZZ.matrix[i].mtype = (clonei.mattype==EXACT) ? IAM_EXACT : IAM_APPROX;
             ZZ.matrix[j].mtype = IAM_MOM;
         }
         else if (clonej.mattype!=PSEUDO && strcmp(clonej.match, clonei.clone)==0) {
             ZZ.matrix[j].high_match = i;
             ZZ.matrix[j].mtype = (clonei.mattype==EXACT) ? IAM_EXACT : IAM_APPROX;
             ZZ.matrix[i].mtype = IAM_MOM;
         }
     }
     if (cnt_cmp >= 1000000) {
         fprintf(stdout, " at clone #%d, compared %.0f, overlaps %d, average %6.3f...\n", 
                  i, (float) total_cmp, cnt_overlap, (float)cnt_overlap/(float) cnt_nodes);
         fflush(stdout);
         cnt_cmp=0;
     }
     
     Zscore_grand2(i);
     
     ZZ.matrix[i].save_grand = ZZ.matrix[i].grand;
return 1;
}


/***************************************************************
        This is the parralelized&optimized full build code
****************************************************************/

/* This does a quick prepass over the matrix - this stuff was removed from
 * the main loop. Sets the first clone to be processed in queueHead.
 */
int Zcreate_node_prepare(int* queueHead)
{
   int i;
   
   struct tmpSz result;
   
   CLONE* clonei;
   
   for (i = 0; i < ZZ.size; i++)
   {
     if (!loadCzfast(&C1z, ZZ.matrix[i].cin))
      {
         ZZ.matrix[i].set = ZNOT_USED;
      }
      else
      {
          ZZ.matrix[i].cb_index = (int *)malloc((C1z.nbands+1) * sizeof(int));

         if (!ZZ.matrix[i].cb_index)
          {
            fprintf(stderr, "Out of memory allocating cb_index.\n");
            return -1;
         }
         
         ZZ.matrix[i].cb_index[0] = ZONE;
         
         clonei = arrp(acedata, ZZ.matrix[i].cin, CLONE);
      

         /* I think this is not used in a complete rebuild.
          * If it is, below comment block will need to be integrated;
          *   this might not be as simple as uncommenting it.
          */
         assert(clonei->mattype == PSEUDO || 
            clonei->match[0] == 0 ||
            clonei->match[0] == ' ');
         
         /*
         if (clonei.mattype != PSEUDO && strcmp(clonei.match, clonej.clone) == 0)
         {
            ZZ.matrix[i].high_match = j;
            ZZ.matrix[i].mtype = (clonei.mattype==EXACT) ? IAM_EXACT : IAM_APPROX;
            ZZ.matrix[j].mtype = IAM_MOM;
         }
         else if (clonej.mattype != PSEUDO && strcmp(clonej.match, clonei.clone) == 0)
         {
            ZZ.matrix[j].high_match = i;
            ZZ.matrix[j].mtype = (clonei.mattype==EXACT) ? IAM_EXACT : IAM_APPROX;
            ZZ.matrix[i].mtype = IAM_MOM;
         }
         */
      }

      if (ZZ.matrix[i].set != ZNOT_USED)
       {
                *queueHead = i;
               queueHead = &ZZ.matrix[i].nextQueued;
      }
   }

   *queueHead = -1;
   
   return 0;
}

/* This loads data from a socket to create the matrix node */
int Zcreate_node_from_socket(int socket, int i, int nodeCount)
{
   int       node;
   u_int32_t j;
   u_int32_t exponent;
   u_int32_t match;
   
   ZNODE* node1;
   ZNODE* node2;

   for (node = 0; node < nodeCount; node++)
   {
      readLoop(socket, &j,        sizeof(j));
      readLoop(socket, &exponent, sizeof(exponent));
      readLoop(socket, &match,    sizeof(match));

      j        = ntohl(j);
      exponent = ntohl(exponent);
      match    = ntohl(match);
      
      Zscore_grand(i, exponent);
      node1 = Znew_node(i, j);
      node1->match = match;
      node1->exponent = exponent;
      
      Zscore_grand(j, exponent);
      node2 = Znew_node(j, i);
      node2->match = match;
      node2->exponent = exponent;
   }
   
   return 1;
}

/* Optimized for a full build - works only after a Zcreate_node_prepare */
int Zcreate_node_parallel(int i)
{
   int  j;
   ZNODE *node, *node2;
   int exponent;
   CLONE clonei, clonej;

   struct tmpCz clone1z, clone2z;
   struct tmpSz result;

   /* It should be impossible to end up doing a buried or unused clone b/c
    * in the work queue they should have been omitted.
    */
   assert(ZZ.matrix[i].set != ZNOT_USED);
   
   /* could be cemap yac with no bands */
   loadCzfast(&clone1z, ZZ.matrix[i].cin); /* always succeeds - prepared */
   
   /* Removed from old loadCzfast, replaced here */
   result.match = result.olap = 0;
   result.prob = 0;

   clonei = arr(acedata, ZZ.matrix[i].cin, CLONE);
   cnt_nodes++;
   
   for (j = i + 1; j < ZZ.size;  j++)
   {
      if (ZZ.matrix[j].set == ZNOT_USED) continue;
      
      clonej = arr(acedata, ZZ.matrix[j].cin, CLONE);
      
      if (LastCB == Select && clonej.selected == 0) continue;
      
      loadCzfast(&clone2z, ZZ.matrix[j].cin);
      /* always succeeds - prepared */
      
      Zsulston_reentrant(&clone1z, &clone2z, &result);
      exponent = ZscoreCpM_reentrant(&clone1z, &clone2z, &result);
      
      if (exponent == 0)
      { /* could be buried with poor overlap */
         if (ZZ.matrix[i].mtype == IAM_CHILD)
         {
            if (strcmp(clonei.match, clonej.clone) == 0)
               exponent = 1;
         }
         else if (ZZ.matrix[j].mtype == IAM_CHILD)
         {
            if (strcmp(clonej.match, clonei.clone) == 0)
               exponent = 1;
         }
      }
      
      if (exponent == 0) continue;

      Zscore_grand(i, exponent);
      node = Znew_node(i, j);
      node->match = result.match;
      node->exponent = exponent;
      
      Zscore_grand(j, exponent);
      node2 = Znew_node(j, i);
      node2->match = result.match;
      node2->exponent = exponent;
   }
   return 1;
}

